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The status of an effort to increase the efficiency 
of calculating transient temperature fields in 
complex aerospace vehicle structures is described. 

The advantages and disadvantages of explicit and 
implicit algorithms are discussed. A promising set 
of implicit algorithms with variable time steps, 
known as the GEAR package is described. Four test 
problems, used for evaluating and comparing various 
algorithms, have been selected and finite-element 
models of the configurations are described. These 
problems include a Space Shuttle frame component, 
an Insulated cylinder, a metallic panel for a thermal 
protection system, and a model of the Space Shuttle 
Orbiter wing. Results generally indicate a 
preference for implicit over explicit algorithms 
for solution of transient structural heat transfer 
problems when the governing equations are "stiff" 
Uypical of many practical problems such as insulated 
metal structures). 


NOMENCLATURE 
C capacitance matrix 
DT time step size 
h^ n-th time step 
K conductivity matri x 

Q thermal load vector 

R residual of the system of equations generated by 
the implicit method 
t time 

t^ n-th time point 
T vector of temperatures 

a thermal diffusivity or coefficient in Adams- 
Moulton formula 

3 coefficient in backward difference method 
INTRODUCTION 

An effort is in progress at the NASA Langley 
Research Center to improve capability to predict and 
optimize the thermal -structural behavior of aerospace 
vehicle structures. The focus of this activity is on 


space transportation vehicles such as the Space 
Shuttle Orbiter. A principal task is to reduce the 
computing effort for obtaining transient temperature 
fields. Current activity is focused on evaluation 
and comparison of explicit and implicit solution 
algorithms. 

In reviewing current literature, a preference 
is evident among researchers for implicit algorithms 
for solution of stiff! sets of ordinary differential 
equations (2-7^). Many engineering analysts, however, 
prefer to use the longer-established explicit algo- 
rithms. A partial explanation for this dichotomy is 
that the full power of the implicit approach has not 
been transferred from researchers to engineering 
analysts. In the explicit algorithms, the time step 
is limited (often severely) in order for the tech- 
nique to be stable. In the implicit algorithms, 
there is no stability-imposed limitation on step size. 
The step size is limited by solution accuracy only, 
so that implicit algorithms can, in general, use much 
larger time steps than explicit algorithms. Because 
a single explicit time step is computationally faster 
than a single implicit time step, the key to the 
advantageous use of implicit algorithms is to use the 
largest possible time step size. As presently 1mple-» 
mented in production thermal analysis computer pro- 
grams, implicit algorithms generally require a user- 
specified fixed time step (8-11 ) . The step size must 
be determined by trial and error. 

The strategy being advocated in the solution of 
large problems by implicit methods is to use algo- 
rithms with variable step size and order and to auto- 
matically select both throughout the solution process 
(12-15). A promising set of algorithms, developed 
for the purpose of implementing the aforementioned 


^ Stiff sets or ordinary differential equations are 
characterized by solutions with widely varying time- 
constants. The typical case is when the solution to 
the homogeneous problem has very small time constants 
compared to those of the forcing function (1_). 
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:strategy, is denoted the GEAR algorithms {13-14). 

A version of the GEAR algorithms well -suited to heat 
transfer analysis denoted GEARIB has been recently 
installed in the SPAR finite-element thermal ana- 
lyzer (8) for testing. 

The purpose of the present paper is to describe 
some ongoing evaluations and demonstrations of the 
use of explicit and implicit algorithms for transient 
thermal analysis of structures using the finite- 
element method. A Shuttle frame test article, an 
insulated cylinder, a metallic multiwall thermal 
protection system panel, and a model of the Shuttle 
Orbiter wing are analyzed using SPAR. Comparisons 
between implicit and explicit algorithms are pre- 
sented. The performance of the GEARIB algorithms 
and especially the value of variable step size and 
order is demonstrated. For benchmark checks, the 
cylinder is also analyzed with the MITAS lumped 
parameter program (Tj). It is a characteristic of 
thermal analysis by finite-element and lumped- 
parameter techniques that modeling affects the stiff- 
ness. Since stiffness is one of the key factors in 
the performance of implicit and explicit algorithms, 
the paper contains a study of the effects of 
modeling on the performance of the explicit and 
implicit algorithms. The present work focuses on 
the implicit and explicit algorithms implemented in 
production programs such as SPAR and MITAS. The 
authors do not evaluate but are aware of and hereby 
recognize recent developments underway which are 
still at the research stage. These include, for 
example, the mixed implicit-explicit techniques (17) 
and the use of quasi -Newton methods to solve the 
nonlinear algebraic equations associated with 
implicit algorithms (V8). 

NATURE OF ALGORITHMS USED IN TRANSIENT THERMAL 
ANALYSIS 

A transient heat transfer problem when 
discretized by finite-element, finite-difference, 
or similar techniques, is governed by the following 
system of equations 

Ct = Q(T,t) - K(T,t)T = F(T,t) T(0) given (1) 

where F is generally a nonlinear function. It is 
usually impractical to obtain an analytical solu- 
tion to eq. (1) so that numerical integration 
methods are used. The simplest numerical integra- 
tion technique is the Euler method which uses the 
first two terms in a Taylor series to predict T at 
time tp^^ as 

•T(t„)*h„C-1 F(T(t„),t„) (2) 

Euler's method is an example of an explicit integra- 
tion technique, so-named because T(tp+i) is given 
explicitly in terms of known quantities. Another 
approach is the backward-difference method which is 
an example of an implicit method. In this approach 

'<Vi> ■ 

■T(t„)*h„C-' F(T(t„„), V,) (3) 

Eq. (3) is a system of implicit equations for 

which is generally nonlinear. The 


explicit algorithm is therefore easier to implement 
but must be bounded to avoid numerical instability 
(unbounded propagation of numerical errors during the 
solution). Implicit techniques are generally stable 
and thus can take larger time steps which are deter- 
mined from accuracy considerations. 

Most practical transient thermal analysis 
problems in flight structures have the following 
characteristics which profoundly affect the choice 
of a solution method: 

(1) The thermal response may be divided into 
regions of slowly and rapidly varying temperatures. 
Steep transients accompany initial conditions or 
sudden changes in the heat load. 

(2) The rapidity of variation of the transient 
portion of the temperature history is proportional to 
the quantity l^/a where L is a characteristic 
conduction length and a is thermal diffusitivity. 
During the transient, time steps much smaller than 
L^/a must be taken no matter what type of integration 
technique is used. 

During a period of slowly-varying temperature 
large time steps may be taken by implicit integration 
techniques but explicit techniques must still use 
time steps which are less than Lva. When l^/a values 
for some elements in the structure are small compared 
to the time scale of the slower temperature variation, 
the problem is stiff. It follows that stiff problems 
are usually best solved by implicit methods. The 
effort involved in solving a system such as eq. (3) 
is usually cost-effective if a small number of large 
time steps are used. 

The Euler method and the backward-difference 
method are presented as representatives of a large 
class of explicit and implicit techniques, respectively. 
Higher-order methods (i.e., multistep) typically use 
more previous information to predict the temperature 
at the current time but the stability properties of 
explicit multistep methods are similar to those of 
the Euler method. Most explicit methods are unstable 
for time steps much larger than L2/a. Accordingly, 
thermal analysis computer programs generally select 
the explicit time step automatically based on the 
stability requirement. For implicit methods, the 
analyst is left to select the implicit time step and 
order without a great deal of guideline information 
and usually several trial runs are needed, There is 
an emerging consensus that the approach to take for 
integrating stiff systems of ordinary differential 
equations would be to use implicit methods which 
automatically select the order and the step size 
based on desired accuracy. One package denoted the 
GEARIB algorithms has these features and is discussed 
next. 

THE GEARIB ALGORITHMS 

Several software packages based on the work of 
Gear have been developed for general use (1_3). The 
package most appropriate for application to finite- 
element thermal analysis is denoted GEARIB. This 
package is intended to solve systems of ordinary 
differential equations of the form 

C(T.t) t = F(T,t) (4) 

The package employs two classes of implicit multistep 
methods, Adams-Moulton and backward difference. For 
nonstiff equations, the Adams-Moulton method of order 
one through twelve is used. This method has the 
general form 
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(5) 


-ALUMINUM STRUCTURE 


T(t, 


n+1 


) = T(t„) + h, 


n Z 
i=0 


where q is the order. For stiff equations, the 
backward difference algorithms of orders one through 
five are used. Jhese algorithms have the general 
form 


i=l 

The coefficients and are given in 

(15). The user selects the class of methods (Adams- 
Mouton or backward differences), and as described 
in (j[3) GEARIB automatically selects the appropriate 
time step and the order based on a user-specified 
error tolerance. 

Use of the GEARIB algorithms is illustrated 
using the backward difference option. Applied to 
eq. (4), eq. (6) gives 

i=l 

-^6o I'l 

This system of nonlinear algebraic equations is 
solved by the modified Newton's method. That is 


where 

i|t1 ■ ' • W 

and J = 9F/8T is the Jacobian of the system at a 
previous time point. Methods used in GEARIB for 
computing J are described in (1_3) and (19) . 

DESCRIPTION OF TEST PROBLEMS AND RESULTS^ 

Insulated Shuttle test frame 

A Shuttle Orbiter frame component analyzed and 
tested under transient heating as described in (20) 
is shown in figure 1 and consists of an aluminum 
frame surrounded by insulation. The principal 
purpose of the study of the configuration as 
discussed in (20) was to evaluate the thermal per- 
formance of the insulation during a simulated Shuttle 
flight. A secondary purpose was to evaluate the 
adequacy of thermal analysis techniques applicable 
to the Shuttle. 


^ Additional details of the test problems are 
given in (19). All calculations were performed on the 
Langley CDC Cyber 173 computer. 



Fig.l Finite-element and lumped parameter models of' 
Shuttle frame 

The lumped parameter model from ( 20 ) consists of 
a two-dimensional section of a symmetric half of the 
structure and contains 118 nodes (see figure 1). The 
unknown temperatures are located at the centroids of 
the lumps. The lumped parameter model was converted 
to a finite-element model for analysis using the SPAR 
program (8). The corresponding SPAR finite-element 
model contains 149 grid points located at the ends 
or corners of the elements. The model contains 148 
elements including one-dimensional elements which 
account for conduction in the aluminum structure and 
radiation across the air gap and two-dimensional 
elements which model conduction in the insulation, and 
across the gap. The difference in numbers of elements 
and grid points is due to the different modeling 
approaches of the two methods. 

Minor modifications were made to the finite- 
element model following the conversion. These con- 
sisted of eliminating or consolidating some extremely 
thin or short finite elements in the aluminum structure 
in order to. reduce the stiffness of the equations and 
to increase the allowable time step for the explicit 
solution algorithm. The properties of the aluminum 
structure are functions of temperature and the prop- 
erties of the insulation are functions of temperature 
and pressure. The pressure dependence is treated in 
SPAR as time dependence using the pressure vs. time 
variation from the trajectory data for the simulated 
flight conditions. The applied heating is specified 
by tabulations of temperatures at the outer surface 
of the insulation. 

The temperature history for the frame was com- 
puted using explicit (Euler) and implicit techniques 
(Crank-Nicholson and backward differences) and GEARIB. 
Comparisons of solution times are given in Table I. 

The explicit procedure using a time step of 0.16 s 
required 1723 s of CPU time. This time step was 
controlled by conduction through most of the aluminum 
elements along the center and front of the frame. 
Solution time using the Crank-Nicholson algorithm 
varied from 475 s to 65 s as the time step was varied 
between T.O and 50 s, The solution times for back- 
ward differences were close to those of Crank- 
Nicholson and are not shown. The GEARIB algorithm 
used time steps from 50 to 170 s and the solution 
time was 54 s. As indicated in Table 1(b), there is 
very little loss of accuracy in either the structure 
or insulation temperatures with increased time step 
size. 
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TABLE I.- PERFORMANCE OF VARIOUS ALGORITHMS FOR TRANSIENT THERMAL ANALYSIS OF SHUTTLE FRAME 

(a) Solution Time Comparison 


Explicit 

! Implicit 1 

Euler 

Crank- 

Nicholson 

1 GEARIB 1 

Time Step (s) 

Solution Time (s) 

Time Step (s) 

Solution Time (s) 

Time Step (s) 

Solution Time (s) 

0.16 

1723 

1 

475 

50-170 

54 



10 

249 





'25 

106 





50 

65 




(b) Effect of Time Step on Accuracy of Implicit Algorithms 


Step Size (s) 

Temp, of Node 309** at 1200 s 

Temp, of Node 

49** at 1200 s 

K 

“F 

K 

“F 

1.0 

442.1 

335.7 

477.0 

398.6 

10.0 

442.0 

335.6 

476.9 

398.5 

25.0 

439.4 

331.6 

475.6 

396.0 

50.0 

437.9 

328.3 

474.8 

394.7 

0.16* 

442.1 

335.7 

477.0 

398.6 

50-170*** 

443.1 

337.5 

477.9 

400.3 


♦Explicit Algorithm 
**See figure 1 
***GEARIB 


TEMPERATURE APPLIED OUTER SURFACE TEMPERATURE (NODE 29) 

Op K ANALYSIS (SPAR EXP. (DT = . 16s). IMP, (DT = 50s) 



Multiwall thermal protection system panel 

The next example problem is one which grew out 
of a study of the thermal performance of a titanium 
multiwall thermal protection system (TPS) panel which 
is under study for future use on space transportation 
systems (21). The configuration as depicted in 
figure 3(aT consists of alternating layers of flat 
and dimpled sheets fused at the crests to form a 
sandwich. The representation of a typical dimpled 
sheet is shown in figure 3(b). For the purpose of 
this analysis, it is assumed that the heat load does 
not vary in directions parallel to the plane of the 
panel* This assumption in addition to the regular 
geometry of the structure leads to the modeling 
simplification wherein only a triangular prismatic 
section of the panel needs to be modeled; fig. 3(a). 

The intersection of this prism with a typical dimpled 
layer is indicated by the shaded triangle in fig. 3(b). 


Fig. 2 Temperature history in outer structural 
surface of Shuttle frame (Node 309) 

The accuracy of the solutions by the various 
techniques is further assessed in figure 2 which 
displays temperature histories at a point in the 
outer layer of the aluminum structure corresponding 
to node 309 (see figure 1). The solid line in 
figure 2 represents the applied temperatures at the 
outer surface of the insulation (node 29). The 
dotted line shows temperatures obtained by the SPAR 
analysis. The SPAR temperatures are plotted as a 
single curve since there is little difference 
between the results. The dashed-dot line shows 
analytical results from the lumped parameter 
analysis of (20) which are also in close agreement 
with the SPAR temperatures. The. circular symbols 
represent test data from (20). The closeness of 
all the results indicates tfiat the models are 
adequate to simulate the temperature history in the 
test article. 




(b) REPRESENTATION OF DIMPLED LAYER 



(C) FINITE ELFMENT MODEL 


Fig. 3 Multiwall thermal protection system panel 
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The finite-element model shown in fig. 3(c) 
contains 333 grid points located on nine titanium 
sheets (five horizontal and four inclined). The 
model contains 288 triangular and quadrilateral 
metal conduction elements, 264 solid air conduction 
elements which account for gas conduction between 
the layers and 544 triangular and quadrilateral 
radiation elements which account for radiation heat 
transfer between adjacent horizontal and inclined 
sheets. Thermal properties of titanium and air are 
functions of temperature. Radiation exchange (view) 
factors were computed and supplied to SPAR using 
the TRASYS II computer program (^). 

The temperature history of the panel in response 
to an imposed transient temperature at the outer 
surface of the panel was computed for 3200 s. 

Results were obtained with SPAR using explicit, 
Crank-Ni Choi son, backward difference and GEARIB 
algorithms. Solution-time comparisons are presented 
in Table II. The explicit algorithm required a time 
step of .007 s. This time step was dictated by 
conduction of heat through the short heat paths 
between the vertices of adjacent triangular layers 
and indicates that this is an extremely stiff 
problem. Required solution time for the explicit 
algorithm was estimated to be 98368 s. 

TABLE II.- COMPARISON OF ALGORITHMS FOR TRANSIENT 
THERMAL ANALYSIS OF TITANIUM MULTIWALL TPS 
(3200 s temperature history) 


Explicit 

Impl icit 

Euler 

Crank- 

Nicholson 

GEARIB 

Time 

Step 

(s) 

Solution 
Time (s) 

Time 

Step 

(s) 

Solution 
Time (s) 

Time 

Step 

(s) 

Solution 
Time (s) 

.007 

98368* 

1 

28412** 

1.0- 

2754 





113 




5 

6352 




^Extrapolated value based on 12296 s for 400 s of 
temperature history 

**Extrapolated value based on 8879 s for 1000 s of 
temperature history 

The Crank-Nicholson solution was carried out 
using time steps of 1 and 5 s which led to solution 
times of 28412 s and 6352 s, respectively. Back- 
ward difference was used with the same time steps 
and had the same solution times. GEARIB took time 
steps ranging between 1.0 and 113 seconds and 
required a solution time of 2754 seconds. This 
example shows again advantages of using implicit 
algorithms in general and the GEARIB algorithms in 
particular for thermal analysis of stiff problems. 

A plot of typical temperature histories for a point 
midway through the panel and the primary structure 
is shown in figure 4 along with the applied outer 
surface temperature. The results were obtained by 
the implicit algorithm with a time step of 5 s and 
are identical to results using a time step of 1 s 
and GEARIB. 




-J 2 
32x10^ 


12 16 
TIME, sec 


20 


Fig. 4 Transient temperatures in titanium multiwall 
TPS panel 


Space Shuttle orbiter wing 

The SPAR thermal model of the Shuttle orbiter 
wing (figure 5) consists of a relatively coarse model 
of the structure (327 grid points) augmented by layers 
of insulation attached to the upper and lower surfaces. 
The structure is modeled by rod, triangular and quad- 
rilateral elements (K21, K31 , K41 SPAR elements). 

The insulation on each surface is modeled by six 
layers of one-dimensional conduction elements (K21). 

Use of these elements neglects lateral heat transfer 
in the insulation--a reasonable assumption since the 
temperature gradients through the insulation are at 
least an order of magnitude greater than the lateral 
temperature gradients. The complete model contains 
2289 grid points, 1400 one- and two-dimensional 
elements in the structure and 1962 one-dimensional 
elements in the insulation. Thermal properties of 
the aluminum structure are temperature-dependent; 
thermal properties of the insulation are temperature- 
and time-dependent. 
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J- 
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Fig. 5 Transient temperatures in Shuttle orbiter wing 
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TABLE III.- COMPARISON OF ALGORITHMS FOR TRANSIENT 
THERMAL ANALYSIS OF SPACE SHUTTLE ORBITER WING 
(4500 s temperature history) 


Explicit 

Implicit 

Euler 

Crank 

-Nicholson 

GEARIB 1 

Time 

Solution 

Time 

Solution 

Time 

Sol ution 

Step 

(s) 

Time (s) 

Step 

(s) 

Time (s) 

Step 

(s) 

Time (s) 

10 

2288 

10 

11730 

1 

cn 

ro 

00 

557 


For the purpose of this analysis, the applied 
heating on the wing is represented by a time- 
dependent temperature applied to the external surface 
of the insulation on the under side of the wing. 

The shape of this curve shown as the solid line in 
figure 5 is roughly indicative of atmospheric reentry 
heating. The temperature history of the wing for 
4500 seconds was computed using the explicit, Crank- 
Nicholson, backward difference and GEARIB algorithms. 
Solution time comparisons are shown in Table III 
along with the time steps used to obtain comparable 
accuracy. The explicit algorithm used a time step 
of 10 seconds--in fact stability requirements 
actually permitted a time step of over 100 seconds 
but the step size was dictated by accuracy and the 
need to periodically update temperature-dependent 
material properties and not by stability requirements. 
The large permitted time step is due to the coarse 
modeling of the structure which did not include the 
thin, high-conducting or radiating elements present 
in the previous models. The implicit algorithms 
(Crank-Nicholson and backward difference produced 
the same results) were used with a time step of 10 s 
and required about five times as much computer time 
as the explicit algorithm. The GEARIB algorithms 
performed very well for this problem. By adaptively 
varying the time step from 1.0 second early in the 
temperature history to as large as 528 seconds toward 
the end, GEARIB required only 557 seconds to complete 
the solution. Figure 5 shows the temperature- 
histories of a point on the structure and a point 
in the insulation 1/5 of the distance through the 
insulation of a typical cross section through the 
wing. The explicit, implicit, and GEARIB algorithms 
produced essentially the same results. 


elements (K41) which receive the heat load and 
radiation elements (R41) which radiate to space. 

Model I contains 800 grid points and 650 elements. In 
model II, the solid elements in the structural layer 
are replaced by quadrilateral elements (K41) in which 
temperatures do not vary through the thickness. This 
is generally a good assumption for thin metal structures. 
Model II has an extra layer of solid elements in the 
insulation in order to preserve the number of grid 
points in the model at 800. In model III, the insula- 
tion is modeled with one-dimensional conductors (K21). 
This model neglects lateral heat conduction but as 
mentioned previously in connection with the Shuttle 
wing model, this effect is small for the class of 
insulated flight structures of interest in the present 
work. 



MODEL I MODEL II MODEL HI 

Fig. 6 Finite-element models of insulated cylinder 

Another aspect of the effect of modeling is 
comparison of results from finite-element and lumped- 
parameter models. To investigate this, the MITAS 
lumped parameter computer program (1^) was applied to 
the analysis of the cylinder. The finite-element 
model I was converted to a lumped-parameter model by 
use of the CINGEN program (^) . The resulting 
lumped-parameter model contained 625 nodes as compared 
to 800 grid points in the finite-element model. Recall 
the unknown MITAS temperatures are located only at the 
centroids of each lump. 


EFFECT OF MODELING ON ALGORITHM PERFORMANCE 

This section of the paper describes a study of 
how modeling details can affect the performance of 
transient solution algori thms--especial ly explicit 
algorithms. Also, the influence of alternate ways 
of including the nonlinear effects of temperature- 
dependent material properties is studied. The 
structure chosen for the study is an insulated 
cylindrical shell shown in figure 6. The cylinder 
is 18 m (720 in.) in length and 4.5 m (180 in.) in 
diameter. The aluminum is 0.25 cm (0.1 in.) thick 
and the insulation is 5.0 cm (2.0 in.) thick. The 
outer surface of the insulation is heated over a 
region which consists of one-third the length and 
half the circumference. 

Three finite-element models are used in the 
study. Due to symmetry, only half the cylinder is 
modeled in each case. In model I, solid (K81) 
elements are used exclusively--39 along the cylinder 
length, 4 around the circumference, and 3 through 
the depth (2 elements in the insulation and 1 in 
the structure). The outer surface has quadrilateral 


TABLE IV.- EFFECT OF MODELING ON SOLUTION TIMES FOR 
INSULATED CYLINDER PROBLEM 


V Problem 
N. and 
NModel 

Algori thm^\^ 

SPAR (Ref. 8) 

MITAS 
(Ref. 16) 

Model 

I 

Model 

II 

Model 

III 

1 umped 
parameter 
model 

Explicit 

10107 

1518 

279 

226 

(time step) 

(.06) 

(2.4-10) 

(3.3-10) 

(10) 

Impl ici t* 
(DT=10 s) 

1880 

1920 

536 

320 

GEARIB 

1779 

1707 

266 


(time step) 

(1.0-83) 

(5-106) 

(2-133) 



*Backward differences and Crank-Nicholson 

The first 2000 seconds of the temperature history 
in the cylinder in response to a time-dependent heat 
load were computed in each model. The explicit 
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(Euler) and implicit (backward difference) algorithms 
were used for all models and in addition GEARIB was 
used for the three SPAR models. Solution times are 
shown in Table IV. Model I is extremely stiff as 
evidenced by the small time step of 0.06 seconds 
required for stability of the explicit algorithm. 

The high stiffness is due to the use of K81 elements 
to model the metal layer. In model II, the stiffness 
has been essentially eliminated by replacing the 
3-D elements modeling the metal by 2-D elements. 

In this model, the explicit technique is faster 
than backward difference and GEARIB. In model III, 
due to low stiffness again, the explicit algorithm 
is faster than the implicit but GEARIB is slightly 
faster than the explicit technique. It is observed 
that in models I and II, GEARIB despite using much 
larger time steps was only marginally faster than 
the implicit method. This is due to the different 
ways of handling the temperature-dependent material 
properties. In the explicit and implicit methods, 
the properties are represented as being piecewise 
constant within time intervals specified by the 
user (by the input quantity TI) in SPAR. Material 
properties are evaluated at the beginning of each 
interval and the conductivity and capacitance 
matrices are regenerated at those times. Results 
for models I, II, and III in Table IV were obtained 
using TI = 20 s. In GEARIB, the material properties 
vary continuously and the residual R must be 
evaluated each time an iteration in solving eq. (8) 
is taken. The residual evaluation is much more costly 
in computer time than the regeneration of the 
conductance and capacitance matrices. This extra 
effort is the price paid for higher accuracy. How- 
ever, this burden only shows up in problems which 
utilize solid (K81) elements due to the extreme 
cost of regenerating the matrices for those elements 
(note model III does not contain K81 elements). A 
way to eliminate the burden (for thermally isotropic 
elements) has been identified and is easily 
implemented. The method is to generate the matrices 
only once for unit values of the appropriate property 
and simply scale the matrices by the property when- 
ever it is updated. 

MITAS computation times are shown in the last 
column of Table IV. Because none of the SPAR models 
is equivalent to the MITAS model in terms of the 
number of unknown temperature or nodal connections, 
no direct comparison of MITAS and SPAR solution 
times is appropriate. However, some trends evident 
in Table IV are noted. The MITAS model is not 
particularly stiff as evidenced by the large time 
step used in the explicit solution technique. SPAR 
models II and III which begin to resemble the MITAS 
model in certain respects are also less stiff and 
favor explicit algorithms. It is noted that the 
way MITAS treats temperature-dependent material 
properties is by the scaling method cited above. 


TEMPERATURE 
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Fig. 7 Effects of choice of algorithm and model 

changes on temperature history of insulated 
cylinder. Model I: all 3-D elements. 

Model II: insulation - 3-D, metal - 2-D. 

Model III: insulation - 1-D, metal - 2-D. 

Figure 7 contains comparisons of temperature 
histories of a point in the cylinder. Model II is 
considered to be the best of the models (recall the 
additional insulation elements used) and thus the 
temperatures represented by the dotted line are 
thought to be the most accurate. These results are 
bracketed by results from model I and MITAS (from 
above) and by model III (from below). There are 
negligible differences between temperatures from the 
implicit and explicit solutions for any given model. 
Results from models II and III are different from 
that of model I because of the extra layer of 
elements through the insulation. The MITAS tempera- 
ture history agrees well with that of model I (on 
which the MITAS model is based) except for some 
differences beginning at 1400 s, 

CONCLUDING REMARKS 

This paper discusses the status of an effort to 
obtain increased efficiency in calculating transient 
temperature fields in complex aerospace vehicle 
structures. Explicit solution techniques which 
require minimal computation per time step and 
implicit techniques which permit larger time steps 
because of better stability are reviewed. A 
promising set of implicit solution algorithms having 
variable time steps and order, known as the GEARIB 
package, is described. Four test problems for 
evaluating the algorithms have been selected and 
finite-element models of each one are described. The 
problems include a Shuttle frame component, an 
insulated cylinder, a metallic panel for a thermal 
protection system, and a model of the Space Shuttle 
Orbiter wing. Calculations were carried out using 
the SPAR finite-element program and the MITAS-lumped 
parameter program. Results generally indicate that 
implicit algorithms are more efficient than explicit 
algorithms for solution of transient structural heat 
transfer problems when the governing equations are 
stiff. Stiff equations are typical of many practical 
problems such as insulated metal structures and are 
characterized by widely differing time constants 
and cause explicit methods to take small time steps. 

As evidenced by their excellent performance in 
solving the test problems, the GEARIB algorithms 
offer high potential for providing increased com- 
putational efficiency in the solution of stiff problems. 
Studies were also made of the effect on algorithm 
performance of different models of the same cylinder 
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test problem. These studies revealed that the 
stiffness of the problem is highly sensitive to 
modeling details and that careful modeling can 
reduce the stiffness of the resulting equations to 
the extent that explicit methods may become 
advantageous. 

REFERENCES 

1 Willoughby, R, A., Stiff Differential 
Systems , Plenum Press, New York, 1974. 

2 Wilson, E. L. , and Nickell, R. E., 
"Application of the Finite Element Method to Heat 
Conduction Analysis," J. Nuclear Engineering and 
Design , Vol. 4, 1966, pp. 276-286. 

3 Wilson, E. L. , Bathe, K. J., and Peterson, 
F. E., "Finite Element Analysis of Linear and 
Nonlinear Heat Transfer," J. Nuclear Engineering 
and Design , Vol. 29, 1974, pp. 240-253. 

4 Hughes, T. , "Unconditionally Stable 
Algorithms for Nonlinear Heat Conduction," Comp. 
Meth. Appl . Mech. Eng ., Vol. 10, 1977, pp. 135-139. 

5 Bathe, K. J. and Khoshogoftaar, M. R., 
"Finite Element Formulation and Solution of Non- 
Linear Heat Transfer," J. Nuclear Engineering and 
Design , Vol. 51, 1979, pp. 389-401. 

6 Hogge, M. A., "Integration Operators for 
First-Order Linear Matrix Differential Equations," 

J. Comp. Meth. in Appl. Mech., Vol. 11, 1976, 

pp. 281-294. 

7 Wood, W. L. and Lewis, R. W. , "A Comparison 
of Time Marching Schemes for the Transient Heat 
Conduction Equation," Int. J. Num. Meth. in Eng ., 
Vol. 9, 1975, pp. 679-689. 

8 Marlowe, M. B. , Moore, R. A., and 
Whetstone, W. D. , "SPAR Thermal Analysis Processors 
Reference Manual, System Level 16," NASA CR-159162, 
1979^ 

9 Trent, D. S. and Welty, J. R. , "A Summary 
of Numerical Methods for Solving Transient Heat 
Conduction Problems," Bulletin No. 49, Oct. 1974, 
Oregon State Univ., Engineering Experiment Station. 

10 Krinke, D. C. and Huston, R. L. , "A 
Critical Evaluation of Computer Subroutines for 
Solving Stiff Differential Equations," Report No. 
UC-ES-1 01 578-8, Oct. 15, 1978, Office of Naval 
Research. 

11 Enright, W. H., Hull, T. E., and Lindberg, 
B., "Comparing Numerical Methods for Stiff Systems 
of O.D.E.'s," m, Vol. 15, 1975, pp. 10-48. 

12 Shampine, L. F and Gear, C. W. , "A User's 
View of Solving Stiff Ordinary Differential 
Equations," SIAM Review , Vo. 21, No. 1, Jan. 1979, 
pp. 1-17. 

13 Hindmarsh, A. C. , "A Collection of Software 
for Ordinary Differential Equations," Rept. No. 
UCRL-82091 , Jan. 1979, Lawrence Livermore 
Laboratory. 


14 Gear, C. W., Numerical Initial Value Problem s 
in Ordinary Differentia Equations . Prentice-Hal l7~ 
Englewood Cliffs, New Jersey, 1971. 

15 Byrne, G. D. and Hindmarsh, A. C., "A 
Polyalgorithm for the Numerical Solution of Ordinary 
Differential Equations," ACM Trans, on Num. Software, 
Vol. 1, No. 1, March 1975, pp. 71-96. 

16 Anonymous, "Martin Interactive Thermal 
Analysis System, Version 1.0," MDS-SPLPD-71 -FD238 
(Rev, 3), March 1972, Martin-Marietta Corp., Denver 
Colorado. 

17 Liu, W. K. and Hughes, T., "Implicit-Explicit 

Finite Elements in Transient Analysis: Implementation 

with Numerical Examples," J. Appl. Mech., Vol. 45, 
1978, pp. 375-378. 

18 Bathe, K. J. and Cimento, A, P. , "Some 
Practical Procedures for the Solution of Nonlinear 
Finite Element Equations," J. Comp. Meth. in Appl. 

Mech. and Eng. , Vol. 22, 1980, pp. 59-85 

19 Adelman, H. M. and Haftka, R. T. , "On the 
Performance of Explicit and Implicit Algorithms for 
Transient Thermal Analysis of Structures," NASA TM 
81880, September 1980. 

20 Gallegos, J. J., "Thermal Math Model of FRSI 
Test Article Subjected to Cold Soak and Entry 
Environments," AIAA Paper 78-1627, 1978. 

21 Jackson, L. R. and Dixon, S. C., "A Design 
Assessment of Multiwall, Metallic Stand-Off, and RSI 
Reusable Thermal Protection Systems Including Space 
Shuttle Applications," NASA TM 81780, 1980. 

22 Jensen, C. L. and Goble, R. G. , "Thermal 
Radiation Analysis System (TRASYS II) Users Manual," 
NASA CR-159273-1 , 1979. 

23 Anonymous, "Creation of Lumped Parameter 
Thermal Models by the Use of Finite Elements," 

NASA CR-1 58944, 1978. 


98 



